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We present the results of a quantum Monte Carlo study of the extended s and the d x 2_ y 2 pairing 
correlation functions for the two-dimensional Hubbard model, computed with the constrained-path 
method. For small lattice sizes and weak interactions, we find that the d x 2_ y 2 pairing correlations are 
stronger than the extended s pairing correlations and are positive when the pair separation exceeds 
several lattice constants. As the system size or the interaction strength increases, the magnitude of 
the long-range part of both correlation functions vanishes. 
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Since the discovery of high-temperature superconduc- 
tivity, the two-dimensional Hubbard model has been the 
subject of an unprecedented level of theoretical activ- 
ity to discover whether it can serve as the paradigm 
for this novel and important phenomenon. Particularly 
with respect to magnetic properties [[I], the physics of 
the model qualitatively represents the behavior of the 
real materials. For example, at half-filling the model is 
an anti-ferromagnetic insulator. Upon doping the anti- 
ferromagnetism rapidly becomes strongly suppressed. 
This behavior is observed in the cuprate superconductors. 
A variety of calculations also predict Q that the doped 
model exhibits an attractive interaction between pairs; 
the s and d x 2_ y 2 symmetries of this attraction are con- 
sistent with the likely symmetries of the experimentally 
measured superconducting gap Yet unobserved, how- 
ever, is convincing evidence that the attractive interac- 
tion leads to a ground state with off-diagonal long-range 
order |l|,0]. In this paper we will present results from 
quantum Monte Carlo (QMC) simulations which suggest 
that the long-range extended s and d x 2_ y 2 pairing corre- 
lations in fact vanish in the thermodynamic limit. 

The fundamental difficulty in deciding whether the 
two-dimensional Hubbard model superconducts is the ab- 
sence of an exact solution. Approximate solutions have 
often been uncontrolled, difficult to benchmark, and con- 
flicting. On several key points, computer simulations 
have provided important information. In fact, the possi- 
ble existence of superconductivity in the Hubbard model 
was suggested by the results of a QMC simulation before 
the discovery of the high temperature superconducting 
materials ||. 

Numerical approaches have, however, had their own 
difficulties, typically being limited to small system sizes, 
high temperatures, and selected electron fillings. Quan- 
tum Monte Carlo methods, for example, experience the 



infamous fermion sign problem Q,|5| , which causes an ex- 
ponential growth in the variance of the computed results 
and hence an exponential growth in computer time as 
the lattice size is increased and the temperature is low- 
ered. While QMC simulations have shown indications of 
pairing correlations, uncertainty has remained because of 
their restriction to relatively small lattice sizes and high 
temperatures [|]||. 

Here, using our new constrained path Monte Carlo 
(CPMC) method §, we discuss the behaviors of the ex- 
tended s and d x 2_ y 2 pairing correlations obtained from 
simulations free of such restrictions. Contrary to stan- 
dard algorithms tms new ground-state (T = OK) 
method exhibits algebraic scaling of computer time with 
system size. It eliminates the exponential growth of vari- 
ances by what we call the constrained-path approxima- 
tion. In a variety of benchmarking calculations , the 
CPMC method has yielded accurate estimates of the en- 
ergy as well as other ground-state observables. 

We considered the following familiar form of the two- 
dimensional Hubbard model on a square lattice of N = 
L x L sites with N a ( cr =t, ].) electrons: 

H=-t^2 (4* c 3° + c ]o c i<?) + U Y1 n *T«a- (!) 

(ij)a i 

we took t — 1 and assumed periodic boundary conditions. 
The pairing correlation function we computed is 

P a (l) = {Al(l)A a (0)), (2) 

where a indicates the nature of pairing. The pair field- 
operator at site I is A a (l) = J2s fa{5)[ci^ci + s j.-Cij,c i+(5 T ], 
where 8 is (±1,0) and (0,±1). For extended s pairing, 
fs(°~) — 1- F° r d x 2_ y 2 pairing, fd{S) is 1 when S — (±1, 0) 
and —1 otherwise. 

To facilitate contact with prior simulations, we also ex- 
amined the "vertex contribution" to the correlation func- 
tion H defined by 
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P a (!) = P«(0 - (At (i)A a (0)) (3) 

The second function on the right is shorthand notation 
for the uncorrelated pairing correlation. For each term 
in P a (l) like (ctcj cjcj), it has a term like (ctc-f)(c|cj_). 

Our numerical method is extensively described and 
benchmarked elsewhere 0. Here we only discuss its 
basic approximation. In this method, the ground-state 
wave function \^o) is projected from a known initial wave 
function \&t) by a branching random walk in an over- 
complete space of Slater determinants \<p). In such a 
space, we can write l^o) = J2<f, x(4>)\4>) ■, where x(4>) > 0. 
The random walk produces an ensemble of \(j>), called 
random walkers, which represent |^o) m the sense that 
their distribution is a Monte Carlo sampling of x(</0- 

To completely specify the ground-state wave function, 
only determinants satisfying (^fo\4>) > are needed; 
hence, |\&o) resides in either of two degenerate halves of 
the Slater determinant space, separated by a nodal plane 
Af that is defined by ( v E , o|0) = 0. The sign problem oc- 
curs because walkers can cross Af as their orbitals evolve 
continuously in the random walk. Asymptotically they 
populate the two halves equally, leading to an ensemble 
that has zero overlap with |\&o)- If Af were known, we 
would simply constrain the random walk to one half of 
the space and obtain an exact solution of Schrodinger's 
equation. Without a priori knowledge of Af, we use a 
trial wave function \^t) and require {SS?t\<I>) > 0. The 
random walk again solves Schrodinger's equation in de- 
terminant space, but under an approximate boundary- 
condition. This is what we call the constrained-path ap- 
proximation. 

The ground-state energy computed by the CPMC 
method is an upper bound. The quality of the calcu- 
lation clearly depends on the trial wave function \^t)- 
Since the constraint only involves the overall sign of its 
overlap with any determinant \<p), it seems reasonable to 
expect the results to show insensitivity to Through 
extensive benchmarking on the Hubbard model, we have 
found that simple choices of this function can give very 
good results . In the calculations reported here we took 
\^t) to be a single Slater determinant. For closed-shell 
electron fillings, we used the free-electron (U = 0) wave 
function. For open-shell fillings, we used unrestricted 
Hartree-Fock (uHF) solutions. For the latter, we have 
found that uHF solutions obtained with low U values 
(< 1), i.e., those resembling free-electron wave functions, 
tend to be good choices for \*S?t) for U up to 8. 

As a calibration of our method, we compare in Table 1 
its prediction for the d x 2_ y 2 pairing correlation function 
to that obtained by an exact diagonalization calcula- 
tion of a 4 x 4 lattice with a closed shell filling of 
2Vj = Ni — 5 and with U = 2 and 4. At each location of 
the pairs, we reproduce the exact result within an error 
of 1 % or less. (Here, in order to compare with the ex- 
act diagonalization data we computed (Ad(i)A^(O)) with 



Arf(Z) = cifJ2sfd{S)ci+si-) 

As a further calibration of our method, we show in 
Fig. 1 the long-range portion of the d x i_ y i pairing cor- 
relation Pd(l) as a function of \l\ for a half-filled 8x8 
system at U = 4. At half-filling, the standard auxiliary- 
field quantum Monte Carlo (AFQMC) method || has no 
sign problem and is exact, and the CPMC method can 
be made exact by removing the constrained path condi- 
tion. In this CPMC calculation, however, we deliberately 
kept the constraint and used for \^t) the uHF solution 
of the system at U = 0.5. In the figure, the computed 
Pd(l) from the CPMC simulation is compared with exact 
AFQMC results. Also shown is the result predicted by 
\^>t)- The inset shows, as a function of electron filling 
(Nf + Ni)/N, the relative difference between the ground- 
state energies calculated by CPMC and AFQMC simula- 
tions. The point indicated by the arrow corresponds to 
the CPMC calculation shown in the main graph. At 1%, 
this difference represents the largest systematic error in 
the CPMC calculation of the energy. With the energy as 
a gauge, the CPMC calculation at 1/2-filling would ap- 
pear to be of the poorest quality; yet, we see that it still 
yields an accurate Pd(l)- The magnitude and range of 
these correlations is comparable to those we now discuss 
for the doped cases. 

Figure 2 shows the long-range part of Pd(l) as a func- 
tion of \l\ for a 12 x 12 lattice at U = 2, 4, and 8. Here, 
the electron filling is 0.85, which corresponds to a closed 
shell case with A T = N± = 61. Figure 2a, the U = 2 
case, shows three different evaluations of this correlation 
function. One is the free-electron prediction for the pair- 
ing. Another is based on definition (|^), and the third 
is the vertex contribution to this definition. Figure 2b 
shows the same set of curves for U = 4 while the inset to 
Figure 2b shows the U = 8 results with the vertex con- 
tribution omitted for clarity. These three sets of curves 
show that Pd(l) is smaller at all three vaules of U than 
the non-interacting case. They also show that increasing 
U causes the long-range correlations, including the vertex 
contribution, to vanish. At U = 8, despite the large error 
bars, the correlations are reduced to simply fluctuating 
around zero. We also see that the vertex contribution 
is a fairly flat function of pair separation up to U = 4. 
This flat region is the "plateau" observed in Ref. jll[] for 
calculations up to U = 2. In this work the Hubbard 
model was studied with next near-neighbor hopping and 
the existence of the "plateau" was attributed to its pres- 
ence. Our results show that the "plateau" behavior is 
no less pronounced in the simple Hubbard model. As U 
increases to 8, however, it vanishes as Pd(l) does. 

In Fig. 3, we address the question of what happens to 
these "long-range" correlations if the lattice size is in- 
creased to 16 x 16. Here, for a closed shell case with the 
same electron filling of 0.85 (JVf = N± — 109), we show 
the CPMC results for U — 2 and 4. First, we notice 
that as in Fig. 2 increasing interaction strength eventu- 
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ally causes the correlations to vanish but now they vanish 
by U = 4. The U = 4 case is shown in Fig. 3b; the accu- 
racy is still sufficient to discern the irregular oscillations 
of Pd(l) around zero. In Fig. 3a the vertex contribution 
is again relatively flat, but nearly zero. Compared to 
Fig. 2a, it has decreased with the increase in lattice size. 
At U = 4 it has in fact vanished and it not shown for 
clarity. We note that we have also carried out calcula- 
tions with a second neighbor hopping and did not find 
any qualitatively different behavior. 

A representative result for the extended s pairing cor- 
relation function P s (l) is shown in Fig. 4 for the same 
system as in Fig. 2. The pairing correlation function 
is shown for the whole range of |Z|. Its short-distance 
magnitude is much greater than that in the tail. We 
mention that the d x 2_ y 2-~w&ve pairing correlation shows 
the same general behavior. In both cases, the short-range 
correlation actually increases as U is increased from zero. 
Hence, the often used integrated pairing correlation func- 
tion, or equivalently the k — (0, 0)-component of P a in 
momentum space, is not a good indicator of supercon- 
ductivity. Comparing the inset with Fig. 2b, we see that 
the long-range extended s-wave pairing is at least an or- 
der of magnitude weaker than the d x 2_ y 2 pairing. Indeed 
it is already fluctuating around zero. 

We also studied pairing correlation functions for other 
electron fillings, lattice sizes, and interaction strengths. 
These results re-enforce those represented above and will 
be reported elsewhere. For example, at open-shell fill- 
ings, even though the calculations experience significant 
increases in variances due to a poorer |^t)'s, results 
through 12 x 12 systems do not appear to show any sig- 
nificant changes compared to the closed-shell ones. 

We reiterate that, due to the constrained-path approx- 
imation, the correlation functions computed here are ap- 
proximate. While the systematic error appears small 
when compared with exact diagonalization and exact 
QMC results, we cannot exclude the possibility that as 
the lattice size increases our systematic error increases 
and an underestimation of the correlations develops. By 
the same token we cannot exclude overestimation. How- 
ever, results like those in Table 1 and Fig. 1, plus a vari- 
ety of other benchmarks M , indicate that our systematic 
error is typically small. In fact, it is often orders of mag- 
nitude smaller than the statistical error in simulations 
using the standard AFQMC method. 

With this very small statistical error, we have pulled 
pairing correlations "out of the noise," and have shown 
examples that for a given system size they disappear as 
the interaction strength increases and for a given inter- 
action strength they disappear as the system size in- 
creases. We note that similar behavior exists for the 
non-interacting problem and the half-filled case. 

We have also computed the lattice size, interaction 
strength, and electron filling dependence of the ground 
state energy, electron momentum distribution, and static 



spin-spin correlation function. We will report these re- 
sults elsewhere. 
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FIG. 1. Long-range behavior of the d x 2_ y 2 pairing correla- 
tion function versus distance for a half-filled 8x8 lattice at 
(7 = 4 computed with |* T ) and by the CPMC and AFQMC 
methods. The inset shows the relative difference between the 
CPMC and AFQMC energies as a function of electron filling. 
The error bars are statistical in origin and mainly associated 
with the AFQMC results. 
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FIG. 2. Long-range behavior of the d x 2_ y 2 pairing correla- 
tion function versus distance for 0.85 filled 12 x 12 lattice at 
U = 2, 4, and 8. This behavior is shown for the free-electron 
and CPMC calculations. Also shown is the vertex contribu- 
tion. 
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FIG. 3. Long-range behavior of the d x 2_ y 2 pairing correla- 
tion function versus distance for a 0.85 filled 16 x 16 lattice at 
U = 2 and 4. This behavior is shown for the free-electron and 
CPMC calculations. Also shown is the vertex contribution. 
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FIG. 4. Extended s pairing correlation function versus dis- 
tance for a 12 x 12 lattice at U — 4 with an electron filling of 
0.85. The inset shows long-range behavior of the free-electron 
case compared with the CPMC result. 



TABLE I. Comparison of the CPMC d x 2_ y 2 pairing correlation function with exact diagonalization results |10| as a function 
of pair separation I — (l x ,l y ). The system size is 4 x 4 with Nf = JVj. = 5. In the CPMC calculations the free-electron wave 



function was 


used for |\&t). Statistical 


errors in these calculations are 


in the last digit and are 


indicated in 


parentheses. 


U 


(0,0) 


(1,0) 


(2,0) 


(1,1) 


(2,1) 


(2,2) 


2 CPMC 


2.0672(2) 


0.0924(1) 


-0.1121(1) 


0.1140(1) 


0.0284(1) 


0.1779(2) 


exact 


2.06693 


0.09223 


-0.11187 


0.11381 


0.02840 


0.17793 


4 CPMC 


2.0635(5) 


0.0876(3) 


-0.0941(3) 


0.1006(4) 


0.0246(2) 


0.1532(6) 


exact 


2.06345 


0.08714 


-0.09422 


0.10013 


0.02453 


0.15302 
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